function [corr_XY,mean_X,std_X,mean_Y,std_Y] = Calc_monthly_corr(stationary_prob_Nz,X_Nz,Y_Nz)
% computes corr(X(t),Y(t)) along with intermediate statistics

    mean_X = sum(stationary_prob_Nz(:).*X_Nz(:));
    mean_Y = sum(stationary_prob_Nz(:).*Y_Nz(:));

    std_X = sqrt(sum(stationary_prob_Nz(:).*(X_Nz(:).^2)) - mean_X^2);
    std_Y = sqrt(sum(stationary_prob_Nz(:).*(Y_Nz(:).^2)) - mean_Y^2);

    corr_XY = (sum(stationary_prob_Nz(:).*X_Nz(:).*Y_Nz(:)) - mean_X*mean_Y)/(std_X*std_Y);

end